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ABSTRACT 

This  thesis  develops  the  theory  for  and  presents  a 
digital  computer  program  capable  of  determining  the  natural 
frequencies  of  a  three-dimensional  piping  system  having 
arbitrary  configuration.   The  analysis  uses  the  method  of 
transfer  matrices.   Piping  hangers,  loops,  and  complex 
branches  (branches  emanating  from  branches)  have  not  been 
included  in  the  analysis.   A  distributed  mass  model  is  used 
for  straight  pipe  sections,  while  mass  is  lumped  for  curved 
sections.   Inclusion  of  shear  deflection  and  rotary  inertia 
is  optional. 

Several  piping  configurations  are  analyzed  using  the 
program;  the  results  are  compared  with  analytical  solutions 
or  values  from  the  literature  to  demonstrate  the  accuracy 
and  integrity  of  the  program. 
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I.   INTRODUCTION 

A.   THE  PROBLEM 

From  an  engineering  standpoint  it  is  important  to  be 
able  to  analyze  and  predict  the  vibration  characteristics 
of  structural  elements  and  assemblies.   The  mathematical 
aspects  of  several  approaches  to  such  analysis  have  been 
known  for  a  long  time,  but  the  volume  of  arithmetic  opera- 
tions dictated  a  gross  simplification  of  complicated  systems 
in  order  to  reduce  the  arithmetic  to  manageable  limits. 
The  development  of  the  electronic  digital  computer  within 
the  last  two  decades,  however,  has  provided  a  means  for 
the  practical  and  accurate  analysis  of  quite  complicated 
vibrating  systems  without  great  simplifications  which  cast 
doubt  on  the  integrity  of  results.   In  turn,  this  has  pro- 
vided a  stimulus  for  further  developments  in  theory  which 
take  advantage  of  computer  capabilities. 

Most  methods  of  analysis  of  actual  systems  in  which  both 
mass  and  elastic  compliance  are  continuously  distributed, 
involve  a  discretization  of  one  sort  or  another;  usually 
the  total  mass  of  the  system  is  considered  to  be  lumped  at 
a  finite  number  of  distinct  points,  the  masses  being  con- 
nected by  massless  spring  elements.   It  is  almost  universally 
assumed  that  the  elastic  characteristics  of  the  system  are 
linear  and  that  the  excursions  of  points  in  the  system  from 
their  equilibrium  positions  are  sufficiently  small  that  no 


geometric  nonlinearities  are  involved.   Such  an  idealization 
leads  to  n  second  order  ordinary  differential  equations. 
Provided  vibrations  are  assumed  to  occur  isochronously 
and  no  damping  is  present,  these  equations  reduce  to  a 
system  of  n  algebraic  equations  in  the  n  amplitudes  and 
frequency  squared. 

The  most  convenient  and  expedient  way  to  cast  this 
problem  is  in  matrix  form.   Thus  engineers  using  these 
methods  deal  with  such  quantities  as  the  mass-matrix, 
the  stiffness  matrix,  the  modal  matrix,  etc. 

On  the  other  hand,  for  certain  simple  cases  of  distributed 
mass  and  elasticity,  it  is  not  necessary  to  lump  the  mass 
into  a  finite  number  of  concentrated  masses.   Instead  the 
vibration  problem  can  be  formulated  as  a  partial  differential 
equation  which  can  be  solved  directly.   Recently,  this  point 
of  view  has  been  extended  for  the  analysis  of  more  complicated 
cases  than  can  be  treated  by  the  conventional  use  of  partial 
differential  equations.   Namely,  the  system  being  analyzed 
is  divided  into  a  number  of  sub-systems,  a  satisfactory 
analysis  existing  for  each  one;  then,  by  a  procedure  which 
will  be  described  in  greater  detail  later,  these  sub- solutions 
are  combined  to  obtain  the  solution  applicable  to  the  original, 
complicated  configuration.   This  procedure,  known  as  the 
transfer  matrix  method,  evolved  from  a  method  described  by 
H.  Holzer  fifty  years  ago  for  the  analysis  of  torsional 
vibrations.   Applications  of  Holzer's  viewpoint  to  the  analysis 


of  flexural  systems  was  made  independently  by  M.  A.  Prohl 
and  N.  0.  Myklestad  about  thirty  years  ago.   The  essentials 
of  these  methods  were  abstracted  by  a  number  of  engineers 
and  mathematicians  who  developed  and  formalized  them  into 
what  is  now  known  as  the  transfer  matrix  method.   The  names 
of  S.  Falk,  K.  Marguerre,  E.  C.  Pestel,  and  F.  A.  Leckie, 
among  others,  are  associated  with  this  development.   A  book 
by  Pestel  and  Leckie  [Ref.  9]  represents  the  fullest  and 
most  accessible  treatment  of  the  subject  available  today. 

The  transfer  matrix  method  is  most  naturally  applicable 
to  the  analysis  of  "chain-like"  structures  composed  of 
elements  which  are  adjoined  at  distinct  points  and  for 
each  of  which  a  satisfactory  vibration  analysis  exists. 
The  simplest  application  is  to  the  case  where  such  a  chain 
of  elements  has  but  two  ends.   However,  it  can  be  applied 
to  structures  having  slightly  greater  topological  complexity, 
but  it  simply  is  not  appropriate  for  the  analysis  of  struc- 
tures which  are  highly  reticulated  or  in  which  major  elements 
are  joined  along  lines  or  surfaces  rather  than  at  distinct 
points . 

Piping  systems  are  substantially  chain-like  and  topologi- 
cally  simple.   Accordingly,  it  was  natural  that  the  transfer 
matrix  method  be  applied  to  the  analysis  of  vibration  in 
piping.   In  this  country  the  first  such  effort  which  was 
reported  was  that  done  by  G.  E.  Fink  in  his  Naval  Postgraduate 
School  thesis  of  1964.   At  about  the  same  time  similar 
efforts  were  being  made  in  Japan  under  the  sponsorship  of 
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the  Japan  Society  of  Mechanical  Engineers.   It  is  known 
that  the  Japanese  have  produced  useful  digital  computer 
programs  for  the  analysis  of  piping  using  the  transfer 
matrix  method,  but  it  has  not  been  possible  to  learn  any 
of  the  details  since  these  programs  remain  under  the  pro- 
prietary control  of  the  Japanese  government. 

Fink  developed  a  program  which  has  convenient  generality 
and  gives  good  results,  but  which  is  limited  to  systems 
lying  entirely  in  a  single  plane.   His  program  is  capable  of 
making  separate  determinations  of  "in-plane"  and  "out-of- 
plane"  vibrations,  which,  for  the  case  of  uniplanar  con- 
figurations, are  uncoupled.   Further  work  on  the  theory 
was  done  by  W.  S.  Baird,  Jr.  and  J.  L.  Simmons  in  their 
Naval  Postgraduate  School  theses;  this  work  mainly  related 
to  systems  having  appreciable  topological  complexity.   In 
his  Naval  Postgraduate  School  thesis,  Y.  S.  Kim  returned 
to  the  development  of  Fink's  program,  adding  several  useful 
features,  including  the  incorporation  of  several  alternative 
mathematical  techniques  for  accurate  solution.   Some  of  these 
techniques  reduced  solution  time  below  that  required  in  Fink's 
original  program;  however,  Kim's  program  is  still  limited  to 
uniplanar  systems. 

Accordingly,  the  writer  undertook  to  develop  a  program 
capable  of  treating  truly  three-dimensional  configurations. 
The  remainder  of  this  thesis  presents  the  analytical  back- 
ground of  this  development,  lists  the  program  which  includes 


rules  for  its  operation,  and  gives  evidence  concerning  the 
high  accuracy  of  the  results. 

It  should  be  made  clear  that  no  claims  are  made  concerning 
the  relative  usefulness  of  the  transfer  matrix  method  as 
compared  to  lumped-mass  methods  which  are  presently  in 
widespread  use  for  analysis  of  piping  vibrations.   In 
general,  the  writer  believes  that  the  results  which  are 
obtained  by  the  transfer  matrix  procedure  are  substantially 
more  accurate  than  can  be  obtained  in  most  cases  using  lumped- 
mass  analysis.   However,  the  delicacy  of  the  calculations 
involved  in  the  transfer  matrix  method  seems  to  imply 
that  only  a  limited  number  of  the  lowest  frequencies  (and 
their  associated  modes)  can  be  obtained  using  present  pro- 
gramming methods  and  computer  hardware.   Greater  word- length 
in  the  computer,  or  availability  of  fast  multiple  precision 
capability,  would  permit  obtaining  a  larger  number  of 
frequencies  and  modes  via  the  transfer  matrix  method. 
Briefly,  the  motivation  for  the  development  reported 
herein  has  been  simply  that  the  transfer  matrix  method  is 
clearly  appropriate  for  the  analysis  of  piping  vibrations 
and  that  the  application  should  be  made.   It  seems  likely 
that  in  some  cases  and  for  some  purposes  the  transfer  matrix 
method  may  be  more  convenient  and  economical  while  in  other 
cases  lumped-mass  methods  may  be  preferable. 
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B.   THE  TRANSFER  MATRIX  METHOD 

Specific  details  of  the  application  of  transfer  matrix 
methodology  to  the  analysis  of  piping  vibrations  will  be 
given  in  subsequent  parts  of  this  thesis.   Background 
material  and  a  general  exposition  of  the  method  is  most 
readily  available  in  the  book  by  Pestel  and  Leckie  [Ref.  9] 
and  it  would  be  pointless  to  include  herein  a  paraphrase 
of  such  material.   However,  for  the  reader  who  is  not  familiar 
with  the  ideas  behind  the  transfer  matrix  method,  a  very 
brief  exposition  is  given  in  this  section. 

An  appropriate  "state  vector"  is  defined  which  presents 
force-type  and  deflection- type  information  capable  of  specify- 
int  the  "state"  at  each  of  several  points  in  the  configura- 
tion at  which  sub-elements  are  joined  together.   These 
quantities  specify  the  configuration  and  internal  force 
system  in  the  structure  at  its  extreme  isochronous,  deflected 
position.   For  example,  in  a  simple  torsional  system,  there 
are  only  two  quantities  in  the  state  vector,  the  angle  of 
rotation  from  the  equilibrium  position  and  the  torque 
transmitted  from  one  element  to  the  next.   The  elastic 
compliance  and  the  inertial  characteristics  of  an  element 
permit  obtaining  equations  relating  the  quantities  in  the 
state  vector  at  the  left  end  of  the  element  to  those  in 
the  state  vector  at  the  right  end  of  the  element.   When 
these  relations  are  put  in  matrix  form  it  is  found  that 
the  state  vector  at  the  right  can  be  written  as  the  product 
of  the  state  vector  at  the  left  premultiplied  by  a  square 
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matrix,  called  a  transfer  matrix,  the  elements  of  which 
represent  the  elastic  and  inertial  properties  of  the  element. 

In  applying  the  method,  a  state  vector  is  chosen  for  the 
left  end  of  the  assembly  of  elements  in  such  a  way  as  to 
satisfy  conditions  of  restraint  at  the  left.   Then,  by  suc- 
cessive multiplication  by  transfer  matrices  representing 
individual  elements  encountered  in  order,  progressing  from 
the  left  end  of  the  assembly  toward  the  right,  one  arrives 
at  a  representation  for  the  state  vector  at  the  right  end. 
The  individual  transfer  matrices  are  functions  of  isochronous 
frequency;  accordingly  the  right  end  matrix  is  a  function  of 
frequency.   When  frequency  is  adjusted,  any  particular 
frequency  which  permits  satisfaction  of  restraint  at  the 
right  end  is  a  natural  frequency  of  the  system.   The  several 
state  vectors  at  the  junction  points  give  the  deflected 
configuration,  or  mode  shape,  corresponding  to  this  frequency, 
when  this  value  is  substituted  in  the  expressions  for  the 
state  vectors. 

For  simple  cases  it  is  possible  to  carry  out  the  multi- 
plications in  literal  form,  keeping  the  frequency  as  an 
unknown  parameter.   Then  the  right  end  restraint  conditions 
provide  a  polynomial  equation  in  the  frequency  which  can 
be  solved  by  customary  methods.   However,  in  relatively 
complicated  cases  it  is  not  feasible  to  retain  the  frequency 
as  a  literal  parameter.   In  these  cases  a  definite  numerical 
value  is  assumed  for  frequency  and  this  leads  to  a  numerical 
measure  of  the  failure  to  satisfy  right  end  restraint 
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conditions.   By  assuming  different  values  of  the  frequency 
and  by  constructing  a  curve  of  measure-of -failure  versus 
frequency,  the  correct  frequencies  can  be  obtained. 

The  transfer  matrices  which  describe  the  individual 
elements  may  be  a  theoretically  complete  and  full  repre- 
sentation, or  they  may  themselves  by  approximate  representa- 
tions.  In  the  application  to  piping  systems  it  is  convenient 
to  represent  a  straight  length  of  pipe  by  a  transfer  matrix 
which  is  theoretically  complete  and  exact.   For  systems 
composed  of  only  straight  lengths  of  pipe  there  is  no 
approximation  ivhatsoever  in  the  transfer  matrix  method 
(other  than  that  implied  by  internal  round-off  in  the  com- 
puter) ;  this  may  be  compared  with  the  fact  that  all  lumped- 
mass  procedures  definitely  imply  an  approximation  and  thus 
result  in  some  theoretical  error. 

Curved  pipe  elements  (elbows  and  bends)  could  be  treated 
exactly  since  adequate  theory  exists  for  the  construction 
of  theoretically  exact  transfer  matrices  with  which  to 
represent  them.   However,  there  are  severe  computational 
difficulties  involved  with  such  exact  representation. 
Accordingly,  with  the  awareness  that  the  total  length  of 
curved  pipe  is  but  a  small  fraction  of  the  total  length 
of  both  straight  and  curved  pipe  in  the  average  configura- 
tion, it  was  decided  to  represent  curved  pipe  by  a  succes- 
sion of  massless  curved  sections  having  proper  elastic 
properties,  alternating  with  appropriately  located  point 
masses.   For  each  of  these  idealized  elements,  it  is 
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possible  to  obtain  the  corresponding  transfer  matrices  with 
no  computational  difficulties.   Thus,  elbows  and  bends  are 
represented  in  much  the  same  way  that  is  done  with  the 
lumped-mass  methods  mentioned  previously.   It  is  believed 
that  the  approximation  introduced  by  this  point-mass 
representation  of  a  small  part  of  a  configuration  is  quite 
tolerable  for  practical  engineering  purposes.   However,  it 
is  certainly  possible  to  modify  the  program  reported  herein 
so  as  to  provide  for  exact  representation  for  curved  ele- 
ments.  This  would  not  involve  major  changes  but  only  the 
development  of  a  subroutine  to  replace  those  given  herein 
for  curved  elements. 


14 


II.   PROGRAM  DEVELOPMENT 

A.   THE  STATE  VECTOR,  TRANSFER  MATRIX,  AND  FREQUENCY 
DETERMINANT 

1 .   State  Vector 

The  state  vector  z.  at  point  i  of  an  elastic  system 

is  a  column  vector  whose  components  are  the  generalized 

displacements,  D.,  and  the  generalized  forces,  F.,  at  that 

point.   In  matrix  notation,  we  have 


Zi  =  {F} 


For  a  beam  free  to  move  in  three  dimensions  with 

longitudinal,  torsional  and  flexural  elasticity,  the  state 

vector   takes  the  form 

r 


z . 
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-  deflection  in  x-direction 

-  deflection  in  y-direction 

-  deflection  in  z-direction 

-  rotation  about  x  -  axis 

-  rotation  about  y-axis 

-  rotation  about  z-axis 

-  axial  force  in  x-direction 

-  shear  force  in  y-direction 

-  shear  force  in  z-direction 

-  Torque  about  x-axis 

-  Moment  about  y-axis 

-  Moment  about  z-axis 


The  order  in  which  the  individual  displacements  and  forces 
occur  in  the  state  vector  is  of  no  particular  importance 
so  long  as  consistency  is  maintained.   Program  VIBREL 
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(for  VIBRation  of  ELastic  piping) ,  the  program  developed 
herein,  employs  the  arrangement  given  above. 
2 .   Sign  Convention  and  Coordinate  System 

A  right-handed  cartesian  coordinate  system  is  used 
to  describe  the  local  generalized  forces  and  displacements, 
the  x-axis  coinciding  with  the  centroidal  axis  of  the  pipe 
as  shown  in  Fig.  2.1. 


*  u     4> 


N 


Figure  2.1:   Vectorial  Representation  of  State  Vector 

Displacement  and  Forces 

A  cut  made  at  any  location  of  the  pipe  will  expose  two 
faces;  the  face  with  its  outward  normal  in  the  positive  x-axis 
direction  is  considered  positive.   The  state  vector  parameters 
are  positive  if  their  vector  representations  coincide  with 
the  directions  prescribed  by  the  cartesian  coordinate 
system  at  the  positive  face  as  in  Fig.   2-1.    Rotation  and 
moment  vectors  are  portrayed  according  to  the  right-hand- 
screw  rule  by  a  double  line  with  an  arrowhead. 

3 .   Transfer  Matrices 

The  transfer  matrix  is  defined  as  the  matrix  which 
relates  two  state  vectors  at  successive  selected  points  of 
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division  in  the  elastic  system. 

z  .  =  U  .  z  -  , 
1     11-I 

where  U.  is  the  transfer  matrix.   The  transfer  matrix  must 
be  square,  and  in  the  case  of  a  spatial  system  must  have  12 
rows  and  12  columns.   When  the  transfer  matrix  relates  the 
flexibility  between  the  points  i  and  i-1  of  a  continuous 
system  it  will  hereafter  be  referred  to  as  a  field  matrix,  V. 
When  a  discontinuity  arises  in  either  force  or  displacement, 
as  with  a  point  mass  in  a  lumped  system,  the  transfer  matrix 
relating  the  state  vectors  on  either  side  of  the  discontinu- 
ity is  called  a  point  matrix,  P.   The  term  "transfer  matrix" 
will  henceforth  be  used  with  regard  to  the  matrix,  U,  which 
includes  mass  and  flexibility  in  the  relationship  between 
the  state  vectors  at  two  distinct  points.   A  complete  discus- 
sion of  transfer  matrices,  with  examples  and  derivations,  can 
be  found  in  Ref.  9. 

4.   The  Elimination  Process  and  Frequency  Determinant 

The  process  of  combining  the  system  transfer  matrices 
to  find  natural  frequencies  can  be  demonstrated  by  a  simple 
planar  example  using  the  Myklestad  method  of  replacing  a 
distributed  mass  beam  with  lumped  masses  connected  by  flexi- 
ble elements  (Fig.  2.2).   Assuming  that  the  field  and  point 
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Figure  2.2:   Lumped  Mass  Representation  of 

a  Planar  Beam 
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matrices  for  a  massless  flexible  beam  and  a  point  mass  are 
already  known,  then 


z,  =»  V,  z  ,  z,  =  Pz,  ,  and  z0  =  V0z.. 
1     lo'l      1 '       2     21 


w 


where  z 


V. 

M 


y 


,   and  superscripts  L  and  R  indicate  left  side 


J  l 


or  right  side,  respectively.   Combining  the  above  yields  the 
relationship 

z0  =  V^PVtZ   =  Uz  (2-1) 

2     2   1  o     o  v    J 

where  U  is  the  system  transfer  matrix.   Written  in  full 

this  relationship  is 

r 


2  L 


11   a12   d13   d14 


21   d22   a23   a24 


31   a32   d33   d34 


41   a42   d43   d44 


w 

V. 
M 


0 


the  a.   being  vibration  functions  of  mass,  stiffness,  and 
frequency.   Application  of  the  boundary  conditions  w   =  w~ 
\\)      =   ty~   =   0  reduces  the  matrix  equation  to  the  following 


a13   a14 


a23   a24 


■ 
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*         ■ 

(2-2) 


Because  the  elements  of  this  square  matrix  are  functions  of 
frequency,  its  determinant  is  called  the  frequency  determi- 
nant of  the  system.   Note  that  this  2X2  matrix  is  a 
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submatrix  of  the  4X4  matrix  appearing  above.   A  non- 
trivial  solution  of  the  equations  2-2  requires  that  the  value 
of  the  frequency  determinant  (A)  equal  zero.   A  value  of  the 
frequency,  a),  which  causes  A  to  vanish  is  called  a  natural 
frequency  of  the  system. 

The  natural  frequencies  could  be  solved  for  directly 
provided  the  expressions  for  the  a.  .  were  kept  in  terms  of 
the  circular  frequency,  oo ,  but  in  practice  these  expressions 
are  forms  much  too  unwieldy  to  be  handled  explicitly. 
Program  VIBREL  evaluates  the  frequency  determinant  of 
order  6X6  for  successive  numerical  values  of  frequency, 
zero  values  of  the  frequency  determinant  corresponding  to 
natural  frequencies  of  the  system. 

Let  the  state  vector  be  composed  of  2r  quantities.   Then 
to  compute  U  in  Eq.  2-1,  two  multiplications  of  2r  X  2r 
matrices  are  required.   Noticing  that  application  of  the 
boundary  conditions  to  the  state  vector  at  point  0  will 
cause  r  of  its  elements  to  be  zero,  we  find  that  r  columns 
of  U  are  multiplied  by  zero,  and  thus  play  no  part  in  the 
calculation.   Kim  [Ref.  6]  introduced  the  concept  of  the 
"state  matrix,"  S,  such  that  the  left  end  state  vector,  zq, 
of  2r  quantities  can  be  written  in  terms  of  a  compressed 
state  vector,  z  * ,  having  only  r  quantities,  viz., 

z   =  Sz  * 
o     o 

where  S  is  the  state  matrix  at  the  left  end,  having  r 

columns  and  2r  rows,  its  elements  consisting  only  of  zeros 
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and  ones.   All  elements  of  the  state  matrix  are  zero  except 
that  in  each  column,  unity  appears  in  the  same  row  position 
as  the  corresponding  non-zero  term  in  the  state  vector, 
boundary  conditions  having  been  applied.   If  the  state 
matrix  is  used  in  the  multiplication  scheme,   then  we  can 
multiply  2r  X  r  matrices  by  2r  X  2r  matrices  and  reduce 
the  calculation  considerably.   The  resulting  2r  X  r  system 
transfer  matrix  requires  only  the  application  of  right  end 
boundary  conditions  to  yield  the  frequency  condition. 

To  illustrate  this  state  matrix  concept  most  simply, 
consider  again  a  two-dimensional  case  in  which  the  state 
vector  consists  of  four  quantities.   The  three  dimensional 
case  is  an  obvious  generalization  of  the  following  discus- 
sion.  Taking  strictly  the  state  vector  approach  we  can 


Figure  2.3: 


Single  Branch  Planar  System 


express  the  relation  between  the  state  vector  at  2  and 
that  at  0  by 


z9  =  U9P1U1z^ 
2     2  1  1  o 


(2-3) 


where  U?  and  U   are  straight  section  transfer  matrices  and 
the  point  matrix  P.  accounts  for  the  effects  of  the  branch. 
Carrying  out  the  matrix  multiplication  and  exhibiting  the 
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system  transfer  matrix  and  state  vectors  which  consider 
bending  only,  we  have 


M 
z 


» 


ail  a!2  a13  al4 


a21  a22  a23  a24 


a31  a32  a33  a34 


a41  a42  a43  a44 


«■            « 

0 

0 

* 

' 

Al 

A2 

(2-4) 


The  boundary  conditions  have  been  applied  to  the  left  end 
state  vector;  the  quantities  A,  and  A~  represent  the  unknown 
quantities  of  V   and  M   respectively.   Application  of  the 
right  end  boundary  conditions  yields  a  frequency  condition 
of  the  same  form  as  Eq.  2-2. 

Using  the  state  matrix  approach  instead,  Eq.  2-3 
becomes 


z~  =  U9P,U1Sz  * 
2     2  1  1   o 


(2-5) 


where  S  is  the  state  matrix 


"o 

0' 

0 

0 

1 

0 

0 

1_ 

contains  the  non- 


zero quantities  of  the  left  end  state  vector.   After  multi- 
plication of  the  matrices  of  Eq .  2-5,  the  number  of  individual 
multiplications  being  considerably  less  than  before,  the 
following  relation  is  obtained: 


••         - 

w 

"  a 
13 

ai4" 

-      i 

Ai 
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r-          — 

a 
23 

a 
33 

a24 
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-. 

A2 
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a43 

a44 
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- 

(2-6) 
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We  see  from  this  result  that  inclusion  of  the  state 
matrix  compresses  the  left  end  state  vector  to  r  quantities 
while  providing  the  same  useful  information  as  the  matrix 
of  Eq.  2-4,  inasmuch  as  the  same  frequency  condition  is 
obtained  after  the  right  end  boundary  conditions  are 
applied.   Also,  if  the  relative  magnitudes  of  the  unknowns 

A,  and  A~  have  been  determined,  the  mode  deflections  and 
forces  at  a  point  can  be  calculated  simply  by  multiplica- 
tion of  the  compressed  left  end  state  vector,  consisting 
of  r  quantities,  by  the  2r  X  r  system  transfer  matrix  up 
to  that  point.   The  appropriate  mode  frequency  is  used  in 
the  computation.   This  will  be  discussed  in  greater  detail 
in  Chapter  3. 

B.  3-D  PIPING  GEOMETRY 

1 .   Working  Points  and  Piping  Nomenclature 

The  geometry  of  the  piping  system  as  it  is  used  in 
program  VIBREL  is  based  on  working  points.   A  working  point 
of  the  system  occurs  where  cross  section  or  properties 
change,  where  branches  join  the  main  member,  where  the 
piping  changes  direction,  and  at  the  ends  of  the  main  member 
and  branches.   These  points  are  measured  in  cartesian  coordi- 
nates from  an  origin  located  arbitrarily  but  fixed  in  space. 
For  purposes  of  this  discussion  we  shall  distinguish  between 
the  two  types  of  working  points  at  which  piping  direction 
changes.   At  a  corner  point  (Fig.  2.4)  the  piping  direction 
changes  with  no  discernible  radius  of  curvature.   At  a 
bend  point  (Fig.  2.5),  which  is  defined  as  the  intersection 
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of  the  tangents  to  the  two  ends  of  a  curved  pipe,  the  radius 
of  curvature  is  finite. 


CORNER 

Point 


-G 


Figure  2.4:   Piping  Corner     Figure  2.5:   Piping  Bend 

The  main  member  is  designated  as  the  longest  continuous 
run  of  pipe  from  which  all  branches  emanate.   Program 
VIBREL  cannot  handle  branches  emanating  from  branches. 

2 .   Unit  Vectors  and  Coordinate  Transformations 

If  a  set  of  mutually  orthogonal  x,  y,  and  z  unit 
vectors  are  formed  at  the  main  member  and  branch  starting 
points,  and  prior  to  and  following  every  every  directional 
change,  we  can  determine,  by  angular  differences  in  succes- 
sive sets  of  unit  vectors,  the  appropriate  coordinate  trans- 
formation angles.   The  local  coordinate  systems  described 
by  these  unit  vectors  are  formed  with  the  piping  system 
in  its  quiescent  configuration. 

Let  us  take,  for  example,  the  two  sections  of  pipe 
pictured  in  Fig.  2.6  determined  by  the  three  working  points 
E,  F,  G, where  F  is  a  corner  point.   Also  let  the  two  sections 
of  pipe,  EF  and  FG,  be  represented  by  vectors  A  and  B,  the 
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Figure  2.6:   Vector  Representation  of  Directional  Change 

components  of  which  have  been  determined  from  the  coordi- 
nates of  E,  F,  and  G.   Then  the  angle  y   can  be  determined  by 
the  following  relation: 


Y |  =  Cos 


-1 


- 

B 

'    A 

|B| 

|  A  | 

The  sign  will  be  determined  later.   Since  the  vectorial 
representation  of  state  vector  quantities  in  Fig.  2.1 
requires  that  the  x-axis  coincide  with  the  centroidal 
axis  of  the  pipe,  let 


|A| 
For  the  purposes  of  succeding  derivations  we  desire  that 
e   be  in  the  plane  of  A  and  B  and  point  outward  from  the 
directional  change.   Hence 


B  x  A 


and   e   =  e   x  e 
y    z    x 


|B  x  A| 

By  the  right  hand  screw  convention  the  rotation  of  amount 
Y  will  always  be  in  the  negative  e   sense.   The  set  of  unit 
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vectors,  0  ,  U  ,  U   are  found  simply  by  rotating  the  initial 

.A.  V  Lt 

unit  vectors  about  the  local  z-axis  by  the  angle  y . 

Now  consider  the  difference  in  orientation  between 

the  0,0,0  vectors  and  those  prior  to  the  next  directional 

x '   y '      z  r 

change  farther  down  the  pipe,0  ,,  U  ,,  and  0  ,.   These  have  been 

J\.  J  L, 

constructed  in  a  manner  similar  to  e  ,  e  ,  and  e  .   Referring 

x '   y '       z  6 

to  Fig.  2.7,  the  angular  difference  between  the 


-  u*  ,  <V 


Figure  2.7:   Unit  Vector  Rotation  about  x-axis 

two  adjacent  sets  of  unit  vectors  is  merely  the  rotation  a 
about  the  local  x-axis.   Alpha  may  be  determined  by  the 
expression 

I  ot  I  =  Cos"1  (0  -0  ') 

This  rotation  is  in  the  positive  (negative)  0   sense  accord- 

ing  to  whether  the  triple  product  (0  xO  , )  .0  f  is  positive 

y   y     x 

(negative) . 

In  this  way,  the  state  vector  is  always  expressed 
in  the  appropriate  local  coordinate  system.   For  the  simple 
system  of  Fig.  2.8 
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Figure  2.8:   Planar  System  with  Directional  Change 


z,  =  U,  z 
1     1  o 


R    TT   L 

z-.  =  U  zn 

1     y  1 


R 


Z2  =  U2Z1 


which  combining  gives 


z0  =  U-U  u\z 
2     2  y  1  o 


where  U   is  the  transformation  matrix  corresponding  to  a 
rotation  y  about  the  z-axis.   For  a  3-D  system  with  the  state 
vector  defined  as  in  section  2.1  this  matrix  is 


U   =  DIAG  [T  ,  T  ,  T  ,  T  ]     ,     T 
Y  y>   y>   y'   y   >  where  T 


Y 


COS  y   SIN  y  0 

-SIN  y   COS  y  0 

0      0    1 


Similarly,  for  the  x-axis  rotations 


U   =  DIAG  [T  ,  T  ,  T  ,  T  ]  ,  where  T 

a  a '   a'   a '   or  '        a 


10      0 

0   COS  a   SIN  a 
0  -SIN  a   COS  a 


C.   3-D  TRANSFER  MATRICES 

1 .   The  Straight  Pipe  Transfer  Matrix 

In  order  to  maximize  accuracy,  a  distributed  mass 
model  was  decided  upon  for  straight  pieces  of  pipe.   The 
program  uses  the  straight  bar  transfer  matrix  cited  in  the 
appendix  of  Ref.  9,  but  revised  to  conform  to  the  ordering 
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of  the  state  vector  described  herein.   The  revised  matrix, 
which  appears  in  Appendix  B,  combines  the  vibration  func- 
tions for  longitudinal,  torsional,  and  flexural  modes  and 
can  include  or  neglect  rotary  inertia  and  shear  deflection. 
2 .   Development  of  the  Curved  Pipe  Field  Matrix 

Due  to  the  complexity  and  increased  machine  time 
involved  in  consideration  of  a  distributed  mass  model  for 
curved  sections  of  pipe,  a  lumped  mass  treatment  was  chosen. 
This,  coupled  with  the  distributed  mass  approach  for  straight 
pipe,  yields  staisfactory  accuracy  combined  with  minimal 
computer  time. 

Given  the  compliance  of  a  piece  of  massless  curved 

pipe  at  its  center  of  curvature,  we  seek  the  field  matrix 

D. 
V  such  that  zR  =  Vz,  (see  Fig.  2.9)  where  z.  =  {p  }.   The 


pJr— ' 


/   \ 

/     \ 

/       \ 

/       \ 

/         \ 
/  \ 


Yl 


Figure  2.9:   Massless  Curved  Pipe 


compliance  matrix  C   is  defined  by  the  relationship 

D   =  C  F  where  D   is  a  vector  of  generalized  displacements 
p    p  p        p 

observed  at  P  and  F   is  some  vector  of  generalized  forces 

P  .......  ^ 

applied  at-  P.   Point  P  is  assumed  to  be  attached  by  a  rigid'" 
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bar  to  one  end  of  the  pipe  while  the  other  end  of  the  pipe 
remains  fixed.   The  compliance  matrix  for  a  pipe  bend  and 
its  derivation  are  presented  in  a  paper  by  J.  E.  Brock 
[Ref.  2]. 

The  analysis  requires  that  we  find  the  compliance 
matrix  at  point  L  which  is  obtained  by  the  appropriate 
rotation  and  translation  of  coordinates  from  point  P  such 

that 

T    T 
CL  ~  BPL  LPL  CP  LPL  BPL 


where   B 


PL 


PL 


PL 


ZL~ZP 


ZP"ZL      yL"y] 


yP"yL      XL"XP 


xp-xL 


for    translation   and 


JPL 


KPL      ° 


0  K 


PL 


'       KPL 


COS(xp,xL)  COS(xp,yL)  COS(xp,zL) 
COS(yp,xL)  COS(yp,yL)  COS(yp,zL) 
C0S(zp,xL)      COS(zp,yL)      COS(zp,zL) 


for  rotation. 


LPL  =  LLP    and    BPL 


;LP 


The  matrix  equation  for  the  generalized  displacements  at 
point  L  can  now  be  written. 

DL  "  CLFL  +  hi   LLR  DR 
where  the  rotation  and  translation  matrices  use  the  above 
notation.   Solving  for  DR  we  arrive  at 
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T    T 
DR  =  LLR  BLR  DL 


T    T 
LLR  BLR  CL  FL 


(2-7) 


The  equation  for  the  generalized  forces  at  the  right  end 
of  the  pipe  has  the  form 

FR  =  "LRL  BRL  FL   * 
If  equations  2-7  and  2-8  are  combined  in  matrix  form,  we 
have  the  result 


(2-8) 


D 


R 


R 


T   T 
LLRBLR 


i_ 


0 


T   T 
LLR  LR  L 


LRLBRL 


DL 

«( 

_FL_ 

By  the  sign  convention  discussed  in  section  2.1,  point  L 
is  located  at  a  negative  face,  so  after  appropriate  sign 
changes  are  taken  into  account,  the  field  matrix  for  the 
bend  becomes 


V  = 


L  TB  T   i 


T   T 
t   r  r 

LLRDLRUL 


0 


LRLBRL 


This  matrix  in  its  explicit  form,  as  used  in  program 
VIBREL,  appears  in  Appendix  B. 

3 .   Development  of  Curved  Pipe  Point  Matrix 

The  mass  of  a  pipe  bend  or  elbow  is  taken  into 
account  through  the  transfer  matrix  procedure  by  lumping  it 
at  the  center  of  mass  and  including  its  effect  on  the  system 
at  point  0  (Fig.  2.10),  the  left  end  of  the  elbow  section. 
The  resulting  point  matrix  relates  only  the  forces  across 
point  0  since  the  deflections  are  continuous. 
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Figure  2-10:   Curved  Pine  Showing  Location 
of  Center  of  Mass 
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A  force  summation  at  point  0  yields  the  following 
scalar  equations  (the  quantities  appearing  in  this  develop 
ment  are  defined  in  Appendix  A,  subscripts  L  and  R  refer 
to  the  left  and  right  of  point  0,  respectively) , 


NR  =  N.  -  mw  (u  -  y0) 


V  =  V  -mw  (v  +  *v 


VzR  =  VzL  -inw  ^w  "  y(f)  "  x^ 


(2-9) 

(2-10) 

(2-11) 


Similarly,  a  moment  summation  at  the  same  point  produces 


Tn  =  TT  -  I    wcf) 
R     L     xo 


M  n  =  M  T  -  I    w   \b 
yR    yL    yo    r 


M  n  =  M  .  -  I    03' 
zR     zL     zo 


(2-12) 
(2-13) 
(2-14) 


Combining  equations  2-9  through  2-14 
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and  the  resulting  point  matrix  equation  is: 


D 


R 


R 


0 


I 


r            -\ 

6 

- 

DL 

6 

FL 

An  inertia  tensor  analysis  was  used  in  the  determination 
of  the  mass  moments  of  inertia  about  point  0.   It  can 
readily  be  shown  that  the  inertia  tensor  at  point  0  is 

!o  ■  ^[ySnf^)ii-SIN2Y(l^li)  +  (f §™2X)3J] 

+  2pAR3[Y-SINY]u  +  pAR3[(-J  +  SIN  |2.)ii 

*(1-C0Sy  -  §™?1HIJ*JI)  -  (|1  -  2SINy  ♦  5IN2Y)33] 

where  it  is  assumed  that  pipe  wall  thickness  is  small 

compared  to  pipe  radius.   (The  symbols  used  here  are  defined 

-?■   =    -}■  -f        =    -t 

in  Appendix  A.)   Then  I    =  i«  I   •  i,  I    =  i  •  I   •  i, 

and  I    =  k  •  I   •  k,  where  i.  i,  and  k  are  unit  vectors 
zo        o     '         '  J  ' 

along  the  local  coordinate  axes  at  point  0.   The  ensuing 
moments  of  inertia  and  point  matrix  are  presented  in  their 
entirety  in  Appendix  B. 

D.   BRANCHED  SYSTEMS 

Any  branch  joining  a  piping  system  has  a  significant 
effect  on  its  dynamic  behavior.   While  the  deflections 
across  the  junction  point  are  continuous,  there  is  a  dis- 
continuity in  the  forces,  the  magnitude  of  which  is 
dependent  on  the  displacement  at  the  junction  point  and 
the  nature  of  the  branch. 
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Consider  the  simple  branch  system  of  Fig.  2.11. 
S   ■■  ^y"  n 

V   J*  B     V___U,„ 


Figure  2.11:   Unit  Vector  Orientation  at  a  Branch  Junction 

Using  the  state  matrix  concept  discussed  in  section  2.1  we 
have  the  folloAving  relationship  between  the  state  vector  at 
point  D  and  the  state  vector  at  point  B  (the  bar  denotes 
branch  coordinate  scheme) . 


DB 
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12  X  1 


R. 


R. 


6X6 


6  X"  6 


-  \6  X  1 


n 


(2-15) 


The  zero  quantities  of  the  state  vector  at  D  have  been 
eliminated  by  application  of  the  boundary  conditions  and 
introduction  of  the  state  matrix,  and  the  branch  transfer 
matrix  has  been  partitioned  into  two  submatrices  P.,  and  P~ 
The  local  coordinate  system  of  the  branch  at  point  B  is 
represented  by  the  unprimed  set  of  unit  vectors.   From  the 
relationship  2-15  it  can  be  seen  that  DR  =  R, z      and 
FR  =  RoZp.,  so  that  when  zn  is  eliminated,  the  equation 
between  the  forces  and  displacements  at  point  B  in  the 
unprimed  coordinate  system  becomes  FR  =  R„R,   DR   (2-16). 
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The  vectors  of  the  generalized  forces  and  displacements 
in  the  branch  system  must  now  be  transferred  to  the  main 
member  system  represented  by  the  unit  vectors  s  in  Fig.  2.11. 
This  can  be  achieved  by  one  z-axis  and  two  x-axis  transfor- 
mations of  the  unprimed  coordinates.   We  first  form  the 
primed  set  of  unit  vectors  in  a  manner  identical  to  that  for 
a  corner  point  (section  2.2),  the  three  working  points 
describing  the  corner  being  E,  B  and  C.   An  x-axis  rotation 
a1  will  align  the  unprimed  with  the  primed  coordinates. 
Next,  a  rotation  y  about  the  z  -axis  will  allow  u  ,  to 
coincide  with  the  centoidal  axis  of  section  BC  resulting  in 
the  double-primed  coordinate  scheme.   Finally,  rotation 
through  the  angle  a  about  the  x"-axis  will  align  the  branch 
system  with  the  main  member  coordinates,  s.   The  total  trans- 
formation of  the  displacements  can  be  expressed  by 


DD=  G   G   G    Dc 

B     a   y   ai   B 


=  GD, 


(2-17) 


where  G   =  DIAG^ 
a 


10     0 
0   COSa   SINa 
0  -SINa   COSa 


10     0 
0   COSa   SINa 
0  -SINa   COSa 


-  » 


similarly  for  G   ,  and 

V 


G   =  DIAG 


COSa   SINa   0 

-SINa  COSa   0 

0      0    1 


COSa   SINa   0 

-SINa  COSa   0 

0      0    1 


Likewise  for  the  forces,  FR  =  GFg  (2-18).   Rearranging 
relations  2-17  and  2-18  and  substituting  in  2-16  yields 
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Noting  that  G 
equation 


-1 


G'\   -  R2Ri1G"1DB   . 


G  we  now  form  the  force  displacement 


-IT   T  T 
FD  =  GR0Rn1G1.1G1G1DT5 
B      2  1   al  y  oi  B 


Referring  to  the  free  body  diagram  in  Fig.  2.12  we  find 
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Figure  2.12:   Generalized  Forces  at  a  Branch  Junction 


that  the  main  member  force  equation  at  point  B  can  be  written 
FR  =  FL  +  FR 

T  _  1   T    T  T 

=  F   +  GLR/G  ,G  G  D„ 

2  1   al  y  06  B  . 


Thus,  since  DR 


L     R 
D   =  D  ,  the  point  matrix  relating  the  state 


vectors  to  the  left  and  right  of  a  branch  junction  point  of 
a  3-D  system  is  given  by 
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The  foregoing  development  assumes  that  there  is  no 
curvature  of  the  main  member  or  branch  at  point  B  and  that 
no  abrupt  change  of  direction  of  the  main  member  occurs  at 
that  immediate  location. 
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For  two  or  more  branches  joining  the  main  member  at 
a  single  point,  the  system  transfer  matrix  just  prior  to 
the  branch  point  is  simply  multiplied  in  turn  by  the  point 
matrix  for  each  of  the  branches. 

E.   COMPUTER  IMPLEMENTATION  OF  TRANSFER  MATRICES 

It  would  be  an  extremely  tedious  job  to  attempt  explicit 
assemblage  of  the  elements  of  the  system  transfer  matrix  in 
terms  of  the  circular  frequency.   Subsequent  solution  for  the 
roots  of  the  expanded  frequency  determinant  would  also  be 
impractical  for  all  but  the  simplest  of  systems.   It  is 
for  these  reasons  that  digital  computation  becomes  a 
necessity . 

Program  VIBREL,  which  incorporates  the  theory  of 
preceding  sections,  was  developed  to  utilize  the  speed  with 
which  the  complete  system  transfer  matrix  can  be  formed  time 
after  time  for  various  values  of  frequency.   A  typical  graph 
of  frequency  determinant  A(w)  versus  circular  frequency  oo 
is  shown  in  Fig.  2.13.   The  values  of  w  at  which  the  graph 
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Figure  2.13:   Graph  of  Frequency  Determinant  vs.  Frequency 
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crosses  the  horizontal  axis,  namely  those  values  of  go  for 
which  A  (go)  =  0,  are  the  natural  frequencies  of  the  system. 

The  program,  presented  in  Appendices  C  through  E,  was 
coded  using  Fortran  language  [Ref.  7]  in  double  precision 
to  minimize  the  effects  of  round-off  and  truncation  errors. 
Since  the  prime  objective  was  to  produce  a  working  program 
of  demonstrable  accuracy  and  reliability,  limited  attention 
was  focused  on  the  niceties  of  programming  intended  to 
reduce  execution  time.   However,  Kim's  state  matrix  concept 
is  included  and  this  effects  a  substantial  saving  of  time 
as  compared  to  other  procedures.   There  is  no  reason  to 
believe  that  the  actual  programming  is  particularly  ineffi- 
cient; however,  it  is  not  unlikely  that  close  scrutiny  of 
program  details  might  point  out  some  places  where  slight 
revisions  would  effect  economies. 

The  transfer  matrices  used  in  computation  of  system 
frequencies  are  those  presented  in  Appendix  B. 

Input  formats  were  designed  to  provide  satisfactory 
simplicity  for  the  user  in  preparing  data  decks,  while 
output  formats  were  chosen  to  check  for  correct  input  as  well 
as  to  provide  a  detailed  summary  of  results.   The  I/O 
formats  are  discussed  more  thoroughly  in  Appendix  E. 
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III.   DISCUSSION 

A.   SCOPE  OF  THE  SOLUTION 

In  its  present  form,  program  VIBREL  is  capable  of 
analyzing  a  piping  system  having  few  restrictions  regarding 
topology.   The  system  must,  however,  be  composed  of  a 
main  member  from  which  can  emanate  as  many  as  fifty  branches. 
The  program  as  listed  in  Appendix  E  can  accommodate  no  more 
than  two  branches  joining  the  main  member  at  a  particular 
point.   In  addition,  there  can  be  no  branches  emanating 
from  other  branches,  nor  can  there  by  any  curvature  of  the 
main  member  or  branch  at  a  branch  junction  point. 

The  total  number  of  sections  capable  of  being  analyzed, 
presently  set  at  one  hundred,  can  be  increased  by  follow- 
ing the  procedure  described  in  Appendix  C.   Likewise  the 
number  of  branches  can  be  augmented  from  its  present  limit 
of  fifty.   For  each  increase  of  one  section,  eighty  addi- 
tional bytes  of  computer  storage  are  needed  above  the  140,000 
bytes  now  required  for  the  entire  program.   For  most  computer 
installations  this  allows  considerable  leeway  in  the  size 
of  the  piping  system  which  can  be  handled  by  the  program. 

To  provide  an  idea  of  the  computer  time  involved  in 
the  frequency  analysis,  the  twelve- section  piping  system 
of  the  sample  problem  in  Appendix  G  required  about  four 
minutes  for  computation  of  four  modes  and  almost  fifteen 
minutes  for  eighteen  modes.   The  IBM  360  computer  compiles 
the  program  in  thirty-nine  seconds. 
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The  usual  assumptions  of  linearity  as  well  as  homogeneity 
and  isotropy  of  the  structural  material  have  been  incorporated 
in  the  analysis.   All  sections  of  the  system  are  assumed  to 
vibrate  isochronously  with  negligible  damping. 

B.   NATURE  OF  ERRORS  AND  INACCURACIES 
1 .   Round-Off  Errors 

As  piping  system  complexity  increases  the  number 
of  multiplications  required  to  form  the  system  transfer 
matrix  grows  proportionately.   Each  time  a  multiplication 
occurs,  there  is  some  round-off  and  loss  of  accuracy;  this 
is  a  function  of  the  significant  figure  capacity  of  the  computer 

FORTRAN,  in  conjunction  with  the  IBM  360  and  double 

precision  arithmetic,  has  the  capability  of  carrying  numbers 

-  78 
of  16  significant  digits  with  an  exponent  range  of  10     to 

7  8 
10   .   Simpler  systems  can  be  checked  for  accuracy  and  it  is 

evident  that  this  significant  digit  capacity  keeps  round- 
off errors  negligible.   Since  with  more  complex  systems 
other  forms  of  errors  creep  into  the  results,  we  have  no 
way  of  predicting  the  exact  contribution  of  round-off 
errors  to  total  inaccuracy  of  frequency. 

It  is  reasonable  to  assume  that  round-off  plays 
a  more  significant  role  in  the  3-D  piping  system  accuracy  than 
in  a  planar  case  with  equivalent  sections  because  of  the 
greater  number  of  coordinate  transformations  that  are  neces- 
sary for  system  description. 
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2 .  Lumped  Mass  Idealization 

For  sections  of  curved  pipe,  lumped  mass  connected 
by  massless  springs  is  used  as  a  model  for  the  real  system. 
This  idealization  leads  to  inherent  inaccuracy  in  the  analy- 
sis of  systems  containing  bends. 

3 .  Zeros  of  the  Frequency  Determinant 

Once  it  has  been  established  that  the  frequency 
determinant  has  crossed  the  zero  axis,  Program  VIBREL  uses 
Mueller's  method  of  successive  bisection  and  inverse  parabolic 
interpolation  to  determine  the  zero.   Care  should  be  exer- 
cised in  considering  as  significant  only  the  same  number  of 
places  as  in  the  acceptability  criterion  which  was  specified 
for  the  solution.   (See  Appendix  C.) 

4 .  Other  Sources  of  Inaccuracies 

The  imperfect  nature  of  the  real  system  and  inevitable 
errors  in  the  measurement  of  working  points  will  render  the 
computer  solution  only  a  good  engineering  approximation  of 
what  can  be  expected  in  actuality.   In  the  discussion  of 
accuracy  in  succeeding  sections,  inaccuracies  caused  by  real 
system  imperfections  will  not  be  considered. 

C.   ACCURACY  AND  INTEGRITY  OF  SOLUTION 
1 .   Solution  Accuracy 

Confidence  in  the  accuracy  of  the  program  solution 
was  established  by  three  comparison  tests,  the  results  of 
which  were  tabulated  in  Appendix  F. 
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a.  Analytical  Comparison 

For  single  straight  sections  of  pipe  a  closed 
form  solution  was  obtained  from  the  governing  fourth  order 
differential  equation  with  the  appropriate  fixed,  free,  or 
ball  joint  end  conditions.   The  largest  variation  between 
the  natural  frequencies  calculated  analytically  and  those 
computed  by  program  VIBREL  was  0.009%. 

b.  Reference  Value  Comparison 

Where  the  exact  closed  form  solution  was  not 
available,  recourse  to  the  literature  provided  some  compari- 
son frequencies.   For  the  simple   curved  pipe  shown  in 
Appendix  F,  Refs.  5,  6,  and  8  contained  a  range  of  frequency 
values  for  the  first  two  modes.   VIBREL  results  for  the 
same  curved  pipe  fell  within  this  range  in  each  of  the  two 
modes  examined,  with  an  average  variation  from  the  reference 
values  of  about  4% .   Because,  for  a  curved  pipe,  program 
VIBREL  includes  bend  flexibility  factor  and  an  accurate 
value  for  shear  distribution  factor  taken  from  Ref  2,  the 
values  computed  are  considered  to  be  at  least  as  accurate  as 
those  contained  in  the  literature. 

For  a  single  branch  system  composed  of  straight 
pipes  only,  comparison  frequencies  for  three  modes  were 
available  from  Ref.  6.   The  largest  difference  between 
reference  and  VIBREL  in  this  case  was  0.06%. 

c.  Dual  Analysis  Comparison 

Dual  analysis  is  based  on  the  fact  that  the  system 
transfer  matrix  is  directionally  dependent.   One  end  of  the 
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system  main  member  is  designated  as  the  starting  point;  then 
the  component  transfer  matrices  are  multiplied  together  in 
sequence  proceeding  toward  the  other  end.   If  an  error 
were  to  occur  in  the  structure  of  one  of  the  transfer 
matrices  or  in  the  method  of  their  combination  to  form  the 
system  transfer  matrix, the  results  from  starting  at  opposing 
ends  could  differ  significantly. 

An  improperly  constructed  state  matrix  would 
also  be  evident  from  a  dual  analysis. 

The  straight  pipe  systems  and  branched  systems 
showed  no  significant  differences  when  subjected  to  dual 
analysis.   The  maximum  difference  observed  in  eighteen 
modes  of  the  complex  system  dual  analysis  was  2.67%, 
while  the  average  difference  was  0.581. 
2 .   Solution  Integrity 

Solution  integrity  was  established  by  observing  that 
the  mode  frequencies  detected  by  the  program  were,  in  fact, 
those  of  the  system  with  no  omissions.   This  was  accomplished 
by  comparison  with  analytical  and  reference  results.   Dual 
analysis  comparison  was  also  used  in  determining  that  the 
same  modes  had  been  detected  in  either  direction.   No  dis- 
crepancies were  noted  in  any  of  the  systems  checked  on  the 
basis  of  these  comparisons. 

Program  VIBREL  includes  in  its  output  a  graph  of 
frequency  determinant  versus  frequency  in  order  to  check 
whether  the  frequency  increment  is  such  that  a  mode  or 
modes  have  been  skipped  in  the  scanning  process.   In  the 
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first  analysis  of  a  branched  system  which  appears  in 
Appendix  F,  the  trend  of  the  output  plot  showed  that  two 
natural  frequencies  in  close  proximity  had  been  overlooked. 
When  the  program  was  rerun  with  the  starting  frequency  incre- 
ment decreased  by  a  factor  of  2,  these  two  frequencies  were 
detected. 

The  maximum  number  of  natural  frequencies  which  can 
be  found  using  VIBREL  is  related  to  the  significant  figure 
capacity  of  the  computer  and  the  physical  configuration 
of  the  piping  system.   As  pointed  out  by  Kim  [Ref.  6], 
at  the  higher  frequencies  the  columns  of  the  frequency 
determinant  approach  the  point  of  being  parallel;  hence  the 
numerical  value  of  the  determinant  includes  differences  of 
large  numbers.   Systems  with  few  components  generally 
exhibit  this  parallelism  at  a  lower  frequency  than  the  more 
complex  ones;  however,  there  is  no  way  of  determining  before- 
hand where  numerical  difficulties  will  be  encountered.   When 
parallelism  is  approached,  a  scattering  of  the  values  of 
frequency  determinant  will  be  noticed  in  the  output  plot. 

Although  the  exponent  limit  on  the  IBM  360  computer 
used  for  testing  the  program  is  ^78,  program  VIBREL  has  been 
coded  to  cease  computation  for  a  particular  problem  and  print 
a  message  when  the  exponent  of  frequency  determinant  exceeds 
+  60. 
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D.   CONCLUSIONS  AND  RECOMMENDATIONS 

Based  on  the  results  of  the  accuracy  and  integrity 
analyses  previously  discussed  and  through  comparison  with 
similar  programs  developed  for  planar  cases,  namely  those 
of  Fink  [Ref  5]  and  Kim  [Ref .  6] ,  program  VIBREL  could  be 
employed  in  its  present  form  as  a  working  tool  for  engineers 
in  the  dynamic  analysis  of  spatial  piping  systems.   The 
fact  that  eighteen  mode  frequencies  were  obtained  for  the 
complex  system  exhibited  in  Appendix  F  implies  that  at 
least  the  first  few  and  usually  the  most  useful  frequencies 
can  be  determined  for  practical  piping  configuration  with 
an  accuracy  consistent  with  engineering  design. 

Minor  modifications  for  adapting  the  program  to  a 
particular  situation  can  be  made  using  the  guidelines  of 
Appendix  C.   With  regard  to  major  revisions,  program  VIBREL 
was  developed  with  solution  accuracy  and  integrity  of  first 
importance.   It  is  recommended  for  future  users  that  changes 
in  the  coding  be  made  for  a  minimum  computer  execution  time 
if  that  becomes  a  prime  requisite.   Avenues  of  investigation 
regarding  numerical  difficulties  at  higher  frequencies  are 
also  open  to  the  future  user.   A  similar  analysis  for  the 
planar  case  was  undertaken  by  Kim  [Ref.  6],  who  utilized 
alternatives  to  the  frequency  determinant  approach. 

Although  the  calculation  of  mode  shapes  of  the  piping 
system  has  not  been  included  in  the  program  because  of  time 
limitations,  an  outline  of  how  to  incorporate  this  feature 
is  discussed  in  the  following  paragraphs.   This  modification 
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to  the  program  is  recommended  as  a  further  important  step 
in  the  dynamic  analysis  of  3-D  piping  systems. 

Let  us  consider  again,  for  simplicity,  the  two-dimensional 
system  of  Fig.  2.3  which  is  repeated  below.   We  found  that 
the  total  system  frequency  condition  could  be  expressed  as: 
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Figure  2.3:   Single  Branch  Planar  System 

A  vanishing  determinant  of  the  coefficient  matrix  yields  a 
natural  frequency.   At  this  frequency  there  is  a  non-trivial 
solution  to  Eqs.  3-1  and  A, /A-  can  be  found.   Normalizing 
A,  and  A-  gives  the  modal  values  of  the  non-zero  quantities 

of  the  left  end  state  vector,  A,*  and  A  * .   A  suitable  nor- 

2     2 
malizing  condition  would  be  A,  +   A?    =    1.   With  these  values 

known  we  proceed  along  the  main  member  calculating  the  modal 

deflections  and  forces  at  each  divisional  point.   This  is 

done  simply  by  multiplying  the  normalized  and  compressed 

left  end  state  vector  by  the  system  transfer  matrix  (2r  x  r) 

In  Fig.  2.3  the  modal  forces  and  deflections  to  the  left  of 

point  1  would  be 


45 


w 


V 


M 
z 

v.   ./ 


lb 


b    bi/i 
13    14 


b22   b24 


b33   b34 


b43   b44 


A  * 
Al 


A  * 

A2 


(3-2) 


In  order  to  permit  unequivocal  interpretation  of  program 
output,  it  would  be  advisable  to  incorporate  a  point  number- 
ing scheme  for  each  point  at  which  modal  forces  and  deflec- 
tions are  desired.   In  addition,  these  forces  and  deflections 
will  be  expressed  in  the  local  coordinate  systems  and  a  method 
should  be  incorporated  for  transforming  them  to  the  global 
system. 

When  a  branch  point,  such  as  1  in  Fig.  2.3,  is  encountered, 
the  branch  transfer  matrix  up  to  point  1  is  formed.   Thus 
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where  B,  and  B-  are  the  unknown  quantities  of  the  state 
vector  at  the  remote  end  of  the  branch.   Since  we  know  that 
the  deflections  are  continuous  across  the  junction  point 
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The  quantities  at  1L  are  given  in  terms  of  the  known  values 
A,*  and  A~*  so  that  the  values  of  B,  and  B~  which  are  intrin^ 
sic  to  the  IB  deflections  may  be  solved  for  through  the 
relations  3-4.   The  known  values  of  B,  and  B~  may  now  be 
used  to  find  modal  forces  and  deflections  in  a  manner  simi- 
lar to  that  used  for  the  main  member. 

Just  to  the  right  of  the  junction  point,  1R,  the  modal 
deflections  are  the  same  as  at  the  left;  however,  the  forces 
must  be  calculated  by  adding  the  1L  forces  to  the  IB  forces 
in  the  global  coordinate  system. 
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APPENDIX  A 
NOTATION  AND  NOMENCLATURE 

The  following  matrix  notation  is  used  in  the  text: 
[  ]  -  matrix 

{  }  -  n  x  1  column  vector 

DIAG  {  }  -  square  matrix  having  the  given  elements 
in  its  principal  diagonal  and  zero 

elsewhere 

T 

-  superscript  denoting  the  transpose  of  a  given 

matrix 

The  following  vector  notation  is  used  in  the  text: 

A  -  vector  A 

|A|  -  length  of  vector  A 

A'B  -  dot  product  of  vectors  A  and  B 

A  x  B  -  cross  product  of  vectors  A  and  B 

Table  A.l  first  lists  the  symbols  used  in  either  the 

test  or  Appendix  B,  with  the  corresponding  computer  program 

variable  identifiers.   Following  these  appear  the  program 

variable  names  which  have  no  counterpart  in  the  text. 

TABLE  A.l 

TEXT    PROGRAM 
SYMBOL   VARIABLE  DESCRIPTION 

A       AREA     cross  sectional  area  of  pipe 

A       A        first  of  two  vectors  defining  a  piping 

directional  change 
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TEXT 

PROGRAM 

SYMBOL 

VARIABLE 

B 

B 

B 

C 

D 

~ex 

UX 

E 

EY 

F 

G 

G 

G 

a 

G2 

G 

Y 

G3 

G  i 
al 

Gl 

i 

Tx 

DIX 

! 

o 

h 

i 

5 

j 

DJ 

k 

k 

BFF 

£ 

DL 

L 

m 

DM 

TABLE  A.l  (Continued) 

DESCRIPTION 

translation  matrix 

second  of  two  vectors  defining  a  piping 
directional  change 

compliance  matrix 

vector  of  generalized  displacements 

unit  vector  in  local  x-direction 

Young's  modulus 

vector  of  generalized  forces 

shear  modulus 

branch  coordinate  transformation  matrix 

branch  coordinate  transformation  matrix 

Branch  coordinate  transformation  matrix 

unit  vector  in  local  x-direction 

mass  moment  of  inertia  of  a  curved  pipe 
about  x-axis 

inei;tia  tensor  about  point  0 

6x6  identity  matrix 

radius  of  gyration  of  pipe  cross  section 

unit  vector  in  local  y-direction 

moment  of  inertia  of  pipe  section 

unit  vector  in  local  z-direction 

bend  flexibility  factor 

length  of  a  straight  section  of  pipe 

rotation  matrix 

mass  of  a  section  of  pipe 
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TEXT 
SYMBOL 

PROGRAM 
VARIABLE 

M 

y 

M 
z 

N 

°6 

P 

R 

R 

f 

RBAR 

Ri 

Rl 

R2 

R2 

rG 

RG 

s 

X 

SX 

t 

T 

T 

u 

u 

0 

X 

UX 

x' 

UXP 

x" 

UXDP 

U 

U,  Ul 

V 

V 

vy 

v„ 

TABLE  A.l  (Continued) 

DESCRIPTION 

moment  about  local  y-axis 

moment  about  local  z-axis 

axial  force  in  local  x-direction 

6x6  null  matrix 

point  matrix 

radius  of  curvature  of  pipe  bend 

mean  radius  of  pipe  section 

submatrix  of  branch  transfer  matrix 

submatrix  of  branch  transfer  matrix 

radius  of  center  of  mass  of  a  curved 
section  of  pipe 

main  member  unit  vector  just  prior  to 
branch  junction 

pipe  wall  thickness 

torque  about  local  x-axis 

deflection  in  x-direction 

unit  dyadic 

unit  vector  in  local  x-direction 

unit  vector  in  local  x'-direction 

unit  vector  in  local  x"-direction 

transfer  matrix 

deflection  in  local  y-direction 

field  matrix 

shear  force  in  local  y-direction 

shear  force  in  local  z-direction 
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TABLE  A.l  (Continued) 


TEST    PROGRAM 
SYMBOL   VARIABLE 


DESCRIPTION 


w 

z 

a 

AL 

al 

AL1 

Y 

GAM 

A 

e 

y 

UU 

6 

PR 

K 

SDF 

P 

RHOM 

4> 

* 

W 

W1,W2 

AA 

ARCL 

BL 

CTX12 

CTZ12 

D 

Dl 

D2 


deflection  in  local  z-direction 

state  vector 

x-axis  coordinate  transformation  angle 

x-axis  coordinate  transformation  angle 

z-axis  coordinate  transformation  angle 

frequency  determinant 

angular  deflection  about  z-axis 

mass  per  unit  length  of  pipe 

Poisson's  ratio 

shear  distribution  factor 

mass  density 

angular  deflection  about  x-axis 

angular  deflection  about  y-axis 

circular  frequency 

array  of  section  properties 

arc  length  of  curved  pipe 

distance  from  bend  noint  to  tangent 
point 

12  x  12  x-axis  coordinate  transformation 
matrix 

12  x  12  z-axis  coordinate  transformation 
matrix 

value  of  frequency  determinant 

outside  pipe  diameter 

inside  pipe  diameter 
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TABLE  A.l  (Continued) 

TEXT    PROGRAM 
SYMBOL   VARIABLE  DESCRIPTION 

DD  frequency  determinant  array 

FN  array  of  natural  frequencies 

GO  gravitational  constant 

ISEC  section  identifier  code 

JBC  branch  boundary  condition  code 

K  boundary  condition  sequence  number 

KB  branch  repeat  discriminant 

LBC  left  boundary  condition  code-main  member 

LBR  branch  point  identifier  discriminant 

M        dimension  of  array  of  frequency  determinant 
values 

MBC  right  boundary  condition  code-main  member 

MCS  curved  subsection  override  discriminant 

MSR  shear  deflection/rotary  inertia  discriminant 

MST  straight  section  test  discriminant 

N  section  number 

NB  vector  of  boundary  condition  codes 

NIB      initial  unit  vector  discriminant  for  a 
branch 

NID      point  identifier  code 

NIT      number  of  iterations  required  for  Mueller's 
method 

NIV      initial  unit  vector  discriminant  for  main 
member 

NK       dual  branch  discriminant 

NMO      number  of  modes  requested  in  analysis 
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TABLE  A.l  (Continued) 


TEXT   PROGRAM 
SYMBOL  VARIABLE 

NN 

NPROB 

NSEC 

NSS 

PX 

PY 

PZ 

RHO 

UA 


UB,UCY 


UST 
WINCR 

WN 

X 

Y 


DESCRIPTION 

iteration  section  number 

number  of  systems  to  be  analyzed  by 
program 

section  number 

number  of  curved  pipe  subsections 

X-coordinate  of  a  working  point 

Y-coordinate  of  a  working  point 

Z-coordinate  of  a  working  point 

weight  density 

point  matrix  for  a  curved  pipe 

field  matrix  for  a  curved  pipe 

straight  pipe  transfer  matrix 

frequency  increment 

natural  frequency 

array  of  frequencies  for  output  plot 

array  of  frequency  determinant  values 
for  output  plot 
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APPENDIX  B 
CATALOG  OF  TRANSFER  MATRICES  USED  IN 


THE    PROGRAM 
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APPENDIX  C 
DESCRIPTION  OF  PROGRAM 

1.   General  Remarks 

Program  VIBREL  is  a  double  precision  FORTRAN  language 
digital  computer  program  which  is  capable  of  performing  a 
frequency  analysis  of  a  three-dimensional  piping  system. 
The  program  accomplishes  this  analysis  by  means  of  the 
transfer  matrix  method  in  which  each  section  of  pipe  is  charac- 
terized by  a  matrix  known  as  a  transfer  matrix.   These  matrices 
which  are  functions  of  frequency  are  multiplied  together 
in  succession  to  form  a  transfer  matrix  which  is  characteris- 
tic of  the  entire  system.   To  combine  satisfactory  accuracy 
with  minimum  computer  time,  a  lumped  mass  model  is  used  for 
curved  sections  of  pipe  while  for  straight  pipe,  a  distributed 
mass  approach  is  employed.   Shear  deflection  and  rotary 
inertia  are  optional  in  the  model. 

The  solution  is  obtained  by  incrementing  the  frequency 
and  forming  the  system  transfer  matrix  and  frequency  determi- 
nant with  each  increment.   The  sign  of  the  frequency  determi- 
nant is  used  to  control  the  search  for  and  convergence  to 
system  natural  frequencies.   The  convergence  acceptability 
criterion  may  be  set  arbitrarily  by  the  user  depending 
on  the  accuracy  desired.   (See  section  3.b  of  this  appendix.) 
As  one  would  expect,  the  greater  the  accuracy,  the  more 
computer  time  required. 
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The  system  is  defined  by  working  points  which  are 
read  into  the  program  with  the  intensive  and  extensive 
properties  of  the  section  of  pipe  preceding  it.   The  details 
for  this  procedure  are  contained  in  Appendix  E.   The  inten- 
sive and  extensive  properties  may  not  vary  along  a  single 
section  of  pipe  but  may  vary  from  section  to  section. 
Projective  complexities,  including  branches  emanating  from 
branches,  cannot  be  handled  directly  by  the  program. 

Three  boundary  conditions,  fixed,  free,  or  ball  joint, 
may  be  specified  by  coded  input  for  the  ends  of  the  main 
member  and  the  remote  ends  of  the  branches.   Alternative 
end  conditions  are  discussed  in  section  3.b  of  this  appendix. 

2.   Program  Structure 
a.   Main  Program 

A  simplified  flow  chart  is  shown  in  Fig.  C.l.   This 
discussion  will  pertain  to  that  diagram. 

Data  which  contains  working  points,  intensive  and 
extensive  properties,  and  boundary  condition  codes  is  read 
into  the  program  as  it  is  required  for  computation.   Local 
coordinate  systems  are  set  up  at  the  starting  point,  prior 
to  and  following  each  directional  change,  and  at  branch 
junction  points.   From  the  working  point  locations  and  unit 
vectors  describing  the  local  coordinate  systems,  section 
lengths  and  coordinate  rotation  angles  are  computed. 

For  computational  purposes  piping  bends  are  sub- 
sectioned  based  on  the  length-to-diameter  ratio  in  the  following 
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Figure  C.l:   Simplified  Flow  Diagram  for  Main  Program 
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manner: 

0  <  L/D  <  1  1  SUBSECTION 

1  <  L/D  £  3  2  SUBSECTIONS 
3  <  L/D  £  6           3  SUBSECTIONS 

If  L/D  is  greater  than  six,  the  number  of  subsections  is 
computed  as  twice  the  number  of  mode  frequencies  sought,  up 
to  a  maximum  of  twelve.   The  user  may,  however,  choose  to 
override  the  number  of  subsections  as  determined  by  the 
program  through  the  use  of  the  parameter  MCS  described  in 
Appendix  E. 

As  the  lengths,  subsections,  and  rotation  angles 
are  computed,  they  are  stored,  together  with  intensive  and 
extensive  properties,  in  an  array.   This  array  is  passed 
via  common  storage  (dotted  line  in  Fig.  C.l)  to  function 
subroutine  FRDET  (FRequency  DETerminant)  which  later  forms 
and  evaluates  the  frequency  determinant  for  a  given  value 
of  frequency. 

When  the  geometry  for  the  entire  system  has  been 
computed,  frequency  iteration  begins.   The  starting  frequency 
is  arbitrarily  set  at  0.1  radians  per  second.   To  compute 
the  starting  frequency  increment,  the  program  constructs 
a  synthetic  straight  pipe,  with  a  length  equal  to  the 
total  length  of  the  main  member,  and  having  outside  diameter, 
wall  thickness,  and  intensive  properties  equal  in  magnitude 
to  a  weighted  average  of  the  main  member  sections.   One 
sixth  of  the  fundamental  frequency  of  this  synthetic  pipe 
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with  fixed  ends  is  taken  as  the  starting  frequency 
increment.   The  system  transfer  matrix  and  frequency  deter- 
minant are  formed  by  function  subroutine  FRDET  for  each 
new  frequency.   When  a  sign  change  is  detected  in  the 
frequency  determinant,  an  iteration  process  using  Mueller's 
method  of  successive  bisections  and  inverse  parabolic  inter- 
polation is  used  for  convergence  to  the  natural  frequency. 
The  convergence  acceptability  criterion  is  taken  as  the  start- 
ing frequency  increment  divided  by  ten  thousand. 

After  a  root  has  been  located,  the  search  begins 
anew  with  the  starting  frequency  equal  to  the  value  of  the 
natural  frequency  plus  one  radian  -per  second.   When  the  second 
mode  frequency  has  been  located,  the  frequency  increment  is 
changed  to  one  tenth  the  difference  between  the  previous 
two  natural  frequencies. 

The  search  process  continues  until  the  specified 
number  of  frequencies  has  been  found  or  until  the  significant 
figure  capacity  of  the  computer  has  been  exceeded.   Then  the 
mode  frequencies  are  printed  and  the  pairs  of  frequencies 
and  frequency  determinant  values  used  in  the  search  are 
sorted  and  plotted, 
b.   Subroutines 

The  program  utilizes  sixteen  subroutines  in  the 

process  of  computation.   A  brief  description  follows: 

ANGLE       Computes  the  sign  and  the  angular 

separation  of  two  sets  of  unit  vectors 
which  have  a  common  axis. 
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COORDX 


Forms  an  n  x  n  (n  is  some  multiple  of 
3  and  <_  12)  coordinate  transformation 
matrix  for  a  rotation  about  a  local 
x-axis . 


COORDZ 


CURMAT 


DETER 


DROOT 


INVERT 


MATMUL 


POINT 


PRPLOT 


SORT 


STATEM 


STRMAT 


UVEC 


Forms  an  n  x  n  coordinate  transformation 
matrix  for  a  rotation  about  a  local  z 
axis  . 

Constructs  the  field  matrix  for  a  curved 
section  of  pipe. 

Computes  the  value  of  an  n  x  n  determi- 
nant; used  to  evaluate  the  frequency 
determinant . 

Performs  successive  bisections  and  inverse 
parabolic  interpolation  to  locate  the 
zeros  of  the  frequency  determinant  after 
a  sign  change  has  been  detected. 

Inverts  the  matrix  R,  to  be  used  in 
computing  the  branch  point  matrix. 


Multiplies  two  matrices  together. 
of  the  matrices  must  be  square. 


One 


Constructs  a  point  matrix  for  a  curved 
section  of  pipe. 

Takes  the  array  of  frequencies  and 
frequency  determinant  values  calculated 
during  iteration  and  plots  them  length- 
wise with  the  printer.   Scaling  is  done 
automatically  and  the  length  of  plot 
may  be  changed  to  suit  the  user. 

Sorts  the  array  of  frequencies  and 
frequency  determinant  values  from 
iteration  for  plotting  on  the  printer. 

Constructs  the  state  matrix  for  the 
appropriate  fixed,  free,  or  ball  joint 
boundary  conditions. 

Constructs  the  transfer  matrix  for  a 
straight  section  of  pipe. 

Computes  a  set  of  orthogonal  unit  vectors 
from  two  pipe  section  vectors  which 
represent  a  directional  change  in  the 
piping. 
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WSTART         Computes  the  value  of  the  frequency 

determinant  for  a  previously  constructed 
synthetic  pipe  equal  in  length  to  the  main 
member  and  having  weighted  average  values 
of  the  properties  of  the  main  member 
sections.   This  frequency  value  is  used 
in  selecting  the  starting  frequency 
increment . 

FRDET  For  a  particular  value  of  frequency, 

this  subroutine  assembles  the  state 
matrices,  transfer  matrices  for  each 
section  of  pipe,  and  the  coordinate 
transformation  matrices.   It  then 
multiplies  them  together  in  sequence 
and  calculates  the  frequency  determinant. 

An  interesting  phenomenon  occurs  in  the  formation  of 
the  frequency  determinant  for  a  single  straight  pipe. 
Because  of  the  symmetry  of  a  pipe  in  the  y  and  z  directions, 
the  two  flexural  submatrices  in  the  straight  pipe  transfer 
matrix  are  identical,  which  causes  the  frequency  determinant 
to  maintain  the  same  sign  on  either  side  of  the  flexural 
natural  frequencies.   This  isolated  case  is  handled  in  sub- 
routine STRMAT  through  the  use  of  the  straight  section  test 
discriminant  and  appropriate  exclusion  of  one  of  the  flexural 
submatrices . 

3.   Remarks  on  Instructions  for  Program  Use 
a.   General  Remarks 

Detailed  instructions  for  program  use  are  contained 
in  the  listing  of  Appendix  E.   It  was  intended  for  these 
instructions  to  be  part  of  the  program  in  order  that  any 
user,  given  only  the  program  listing  or  the  computer  card 
deck,  could  employ  the  VIBREL  package  with  no  prior  knowledge 
of  transfer  matrices  and  without  perusal  of  this  thesis. 
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Some  clarification  of  the  instructions  for  program 
use  is  justified  for  cases  where  properties  or  geometry 
change  within  a  piping  bend.   First,  consider  the  simple 
section  of  piping  shown  in  Fig.  C.2.   Assuming  that  no  change 


.T>-  O 


Figure  C.2:   Piping  Bend  with  No  Change  in  Properties 

in  cross  section  or  properties  occurs  between  points  a  and 
d,  the  only  points  which  must  be  listed  as  working  points 
for  program  input  are  a,  o,  and  d.   Using  these  points  and 
the  radius  of  curvature,  program  VIBREL  will  calculate 
lengths  ab  and  cd  as  well  as  the  arc  length  and  angle  of  the 
bend . 

Now  consider  the  section  of  piping  of  Fig.  C.3  where 
the  cross  section  or  properties  change  at  point  c  but  remain 


Figure  C.3:   Piping  Bend  with  Change  of  Properties 
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constant  between  points  a  and  c  and  between  points  c  and  e. 
In  this  instance  points  a,  o, ,  o~,  and  e  must  be  listed  as 
working  points  for  the  program.   VIBREL  calculates  straight 
lengths  ab  and  de  as  well  as  the  information  required  for 
curved  lengths  be  and  cd.   It  should  be  obvious  from  this 
example  that  for  the  general  case  of  two  curved  sections  occur- 
ring sequentially  that  the  intersection  point  (point  c  in 
Fig.  C.3)  need  not  be  listed. 

If  the  geometry  input  is  such  that  negative  lengths 
are  computed,  an  error  message  is  printed.   One  half  inch 
leeway  is  allowed  for  negative  lengths  to  take  into  account 
any  errors  which  may  occur  in  measurement, 
b.   Program  Flexibilities 

(1)  Approximating  Smarl  Piping  Accessories 

Small  valves,  flanges,  couplings  and  other  piping 
accessories  whose  center  of  mass  is  relatively  close  to  the 
centroidal  axis  of  the  pipe  can  be  approximated  using  a 
straight  section  of  pipe.   The  length,  mass  density,  diameter 
and  wall  thickness  can  be  artificially  varied   to  give  the 
mass  of  the  item.   The  resulting  variation  in  stiffness 
between  an  actual  flange  and  artificial  pipe,  for  example, 
usually  will  not  have  a  significant  effect  on  the  results 
because  of  its  short  length. 

(2)  Approximating  Alternative  End  Conditions 

When  the  ideal  boundary  conditions  are  applied 
to  the  state  vector  on  the  end  of  a  system,  it  must  have 
six  zero  and  six  non-zero  elements.   Actual  field  conditions 
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may  dictate  that  this  is  not  the  case  because  of  flexible 
mountings  or  connection  to  machinery.   A  good  approximation 
in  these  circumstances  may  be  made  by  attaching  one  end  of 
an  artificial  straight  pipe  to  the  system  piping  keeping  its 
other  end  fixed.   The  mass  of  the  artificial  pipe  should  be 
small  with  the  extensive  properties  and  elastic  modulus 
chosen  to  approximate  the  stiffness  of  the  mounting. 

The  program  can  then  use  the  ideal  boundary 
conditions  with  little  loss  of  accuracy  from  the  real  situation. 
(3)  Program  Modifications 

In  some  instances,  a  user  may  desire  to  change 
such  quantities  as  frequency  increment  or  input  format  to 
suit  his  particular  needs.   For  this  reason,  Table  C.l 
lists  the  more  important  program  variables  and  what  card  or 
cards  to  make  these  changes  to. 


TABLE  C.l 
Modification 


Card         Variables 
Numbers       Affected 


Maximum  number  of  points  for        2890  x(400),  y(400) 

print  plot 

Maximum  number  of  piping  sections    2910  AA(100,10) 

in  analysis 

Maximum  number  of  branches  in       2920  NB(50) 

analysis 

Input  format  3210-50 

Output  format  (Input  Data)  3390-420 

Output  format  (Properties  and       7750-60 
Geometry) 

Output  format  (Natural  8740 

Frequencies) 
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TABLE  C.l  (Continued) 

Modification 

Starting  frequency  increment 

Starting  frequency 

Convergence  acceptability 
criterion 

Maximum  iterations  to 
convergence  in  DROOT 

Starting  Frequency  after 
location  of  natural  frequency 

Frequency  increment  after 
second  mode 

Number  of  lines  of  print  plot       14730         RY 


Card 

Variables 

Number 

Affected 

8650 

WINCR 

8760 

Wl 

9120 

EPS 

3670 

NIT 

9330 

Wl 

9320 

WINCR 
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APPENDIX  D 
PROGRAM  FLOW  DIAGRAMS 

1.  General  Remarks 

Flow  diagrams  are  included  for  the  main  program  and 
function  subroutine  FRDET.   Subroutines  DETER,  DROOT, 
INVERT,  and  PRPLOT  are  standard  package  subroutines  with 
minor  modifications  made  to  accommodate  them  to  program 
VIBREL.   Since  their  equivalent  or  an  alternative  routine 
may  be  substututed  for  them,  flow  diagrams  are  not  included. 
Adequate  comment  statements  are,  however,  contained  in  these 
subroutines  to  define  their  structure.   Other  subroutines 
for  which  no  flow  diagrams  appear,  incorporate,  at  most,  two 
decisions  in  their  structure*  - 

2.  Main  Program  Flow  Diagram 
Pages  6?  through  76. 

3.  Function  Subroutine  FRDET  Flow  Diagram 
Pages  77  through  79. 
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APPENDIX  F 
PROGRAM  ACCURACY  AND  INTEGRITY  ANALYSIS 

1.   General  Remarks 

This  appendix  tabulates  by  mode  frequency  (radians/second) 
the  results  of  the  accuracy  and  integrity  analyses  of  various 
piping  systems  from  simple  to  complex. 

The  comparison  frequency  values  for  the  simple  straight 
section  systems  were  obtained  by  solution  of  the  differential 
equations  for  longitudinal,  torsional  and  flexural  vibrations 
with  application  of  the  appropriate  boundary  conditions.  These 
solutions  were  evaluated  on  the  digital  computer  with  accuracy 
to  at  least  as  many  places  as  listed  in  the  tabulation.   Evalu- 
ation of  the  comparison  solution  and  VIBREL  solution  were 

2 
carried  out  using  a  value  of  32.174  FT/SEC   for  the  gravita- 
tional constant  and  a  value  of  tt  to  15  significant  figures. 

Comparison  values  for  the  curved  pipe  analysis  were  obtained 
from  three  sources:   Fink  [5],  Kim  [6],  and  Ojalvo  [8].   Those 
extracted  directly  from  [5]  and  [6]  were  compiled  using  a  lumped 
mass  model  and  digital  computation.   From  graphs  and  formulas 
(based  on  the  theory  of  clamped  ring  segments)  appearing  in  [8], 
the  other  natural  frequency  values  for  the  curved  pipe  were 
calculated . 

Several  comparison  frequencies  for  a  single  branch  system 
were  obtained  from  [6], 

Shear  deflection  and  rotary  inertia  were  neglected  except 
where  noted. 
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2.   SYSTEMS  ANALYSIS 

a.   STRAIGHT  PIPE  WITH  VARIOUS  END  CONDITIONS 


END  CONDITIONS 
FIXED-FIXED 


COMPARISON  RESULTS 

MODE 

1ST 
2ND 
3RD 
4TH 
5TH 
6TH 
7TH 
1ST 
8TH 
9TH 
1ST 


FLEXURAL 

FLEXURAL 

FLEXURAL 

FLEXURAL 

FLEXURAL 

FLEXURAL 

FLEXURAL 

TORSIONAL 

FLEXURAL 

FLEXURAL 

LONGITUDINAL 


SYSTEM  PARAMETERS: 

DIAMETER  (IN) 

WALL  THICKNESS  (IN) 

LENGTH  (IN) 

MODULUS  OF  ELASTICITY  (PSI) 

POISSON'S  RATIO       ^ 

WEIGHT  DENSITY  (LB/ FT-3) 


2.0 

.125 

200. 

30  x 

.25 

490. 


10 


ANALYTICAL 


75.10 
207.0 
405.8 
670.9 
1002. 
1399. 
1863. 
2007. 
2393. 
2990. 
3174. 


461420 

2876110 

5914940 

0  57  898  0 

21750380 

79137910 

62757630 

8325054 

7260868 

0869112 

6619386 


VIBREL 


75.1 
207. 
405. 
670. 
1002 
1399 
1863 
2007 
2393 
2990 
3174 


0461590 

02876110 

85914210 

90577600 

.21750410 

.79140510 

.62799710 

.8275305 

.7157896 

.3653337 

.6560265 


%    DIFFERENCE 

.00000 
.00000 
.00000 
.00000 
.00000 
.00000 
.00002 
.00025 
.00043 
.00931 
.00019 


DUAL  ANALYSIS 

RESULTS : 

STARTING  AT 

STARTING  AT 

MODE 

LEFT  END 

RIGHT  END 

%    DIFFERENCE 

1 

75.1046142 

75.1046142 

.00000 

2 

207.0287611 

207.0287611 

.00000 

3 

405.8591494 

405.8591494 

.00000 

4 

670.9057898 

670.9057898 

.00000 

5 

1002.2175038 

1002.2175038 

.00000 

6 

1399.7913791 

1399.7913791 

.00000 

7 

1863.6275763 

1863.6275763 

.00000 
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END  CONDITIONS 
FIXED-FREE 


SYSTEM  PARAMETERS 
SAME  AS  (1) 


COMPARISON  RESULTS: 

MODE  ANALYTICAL 


11.8028657 

73.9673211 

207.1106480 

405.8541957 

670.9060651 

1399.7913799 


VIBREL 

11.8028696 

73.9673092 

207.1106397 

405.8541953 

670.9060651 

1399.7913885 


;  DIFFERENCE 

00003 
00002 
00000 
00000 
00000 
00000 


ANALYSIS  RESULTS 


STARTING  AT 
LEFT  END 

11.8028696 

73.9673092 

207.1106597 

405.8541953 

670.9060651 

1399.7913885 


STARTING  AT 
RIGHT  END 

11  .8028696 

73.9673092 

207.1106397 

405.8541953 

670.9060651 

1399.7913885 


;  DIFFERENCE 

00000 
00000 
00000 
00000 
00000 
00000 
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T2L 


END  CONDITIONS: 
FIXED-PINNED 


SYSTEM  PARAMETERS 
SAME  AS  (1) 


COMPARISON  RESULTS: 

MODE  ANALYTICAL 


51.7571905 

167.7264474 

349.9478448 

598.4315201 

913.1775124 

1003.9162520 


VIBREL 

51  .7571773 

167.7264180 

549.9478460 

598.4315226 

913.17753170 

1003.9164796 


;  DIFFERENCE 

00003 
00002 
00000 
00000 
00000 
00002 


DUAL  ANALYSIS  RESULTS 


MODE 

1 
2 
3 
4 
5 
6 


START INC  AT 
LEFT  END 

51.7571773 

167.7264180 

349.9478460 

598.4315226 

913.17753170 

1003.9164796 


STARTING  AT 
RIGHT  END 

51.7571773 

167.7  26418  0 

349.9478460 

598.4315226 

913.17753170 

1003.9164796 


DIFFERENCE 

00000 
00000 
00000 
00000 
00000 
00000 
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L^«        \~>  KJ  I  \.    \     I-*  LJ  IXf 


(-/ZO   6>,  So  ,o) 


C-3S.3,-^5.3.0)V 


(35.3,-35.3,0) 


END  CONDITIONS 


FIXED- FIXED 


SYSTEM  PARAMETERS: 

PIPE  DIAMETER  (IN) 
WALL  THICKNESS  (IN) 
INCLUDED  ANGLE  OF  ARC  (DEG) 
MODULUS  OF  ELASTICITY  (PSI) 
POISSON'S  RATIO       - 
WEIGHT  DENSITY  (LB/FT-3) 


2.0 

.125 

270 

30  x  106 

.25 

556 


COMPARISON  RESULTS 

a.   MODE  1 

SOURCE 

VIBREL 
OJALVO  [8] 
FINK  [5] 
KIM  [6] 


FREOUENCY 


37 
38 
38 


2234058 
5930983 
8773520 


35.3639740 


VARIATION  WITH 
VIBREL  VALUE 

3.68 

4.44 

-5.00 


b. 


MODE  2 

SOURCE 

VIBREL 
OJALVO 
FINK 
KIM 


FREQUENCY 

72.1914615 
69.9962862 
73.3297307 
66.6993490 


%    VARIATION  WITH 
VIBREL  VALUE 

-3.04 

1.58 

-7.61 
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c.   BRANCHED  SYSTEMS 

(1)  SINGLE  BRANCH  SYSTEM 


(£9.3,70.7,0) 


7(0,0,0) 


(/00,Ot  o)  (200,0,  o)  I 


SYSTEM  PARAMETERS  (ALL  3  SFCTIONS) : 


DIAMETER  (IN) 

WALL  THICKNESS  (IN) 

LENGTH  (IN) 

MODULUS  OF  ELASTICITY  (PSI) 

POISSON'S  RATIO 

WEIGHT  DENSITY  (LB/FT3) 


2.0 

.125 

100. 

30  x 

.25 

556. 


10' 


COMPARISON  RESULTS: 

MODE  KIM  [6] 


194.25470558 
278.61078241 
628.59556355 


VIBREL 

194.2888499 
278.6390155 
628.9417782 


%    DIFFERENCE 

.01758 
.01013 
.05508 


DUAL  ANALYSIS  RESULTS 


MODE 

1 
2 
3 
4 
5 
6 
7 
8 
9 
10 


STARTING  AT 
LEFT  END 


67. 
194 
198 
242 
278 
281 
370 
628 
634 
689 


6456627 

.2888474 

.0817333 

.0773169 

.6390155 

.9864461 

.7949928 

.9417782 

.7478785 

.7373111 


STARTING  AT 
RIGHT  END 

67.6456627 

194.2888474 

198.0817333 

242.0773169 

278.6390155 

281.9864461 

370.7949928 

628.9417782 

634.7478786 

689.7373111 


%    DIFFERENCE 

.00000 
.00000 
.00000 
.00000 
.00000 
.00000 
.00000 
.00000 
.00000 
.00000 
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(2)  DUAL  BRANCH  SYSTEM 
(40,30,0) 


END  CONDITIONS 

MAIN  MEMBER 
BRANCH  #1 
BRANCH  #2 


FIXED-FIXED 

PINNED 

FREE 


<<?o,-3oto) 


SYSTEM  PARAMETERS: 


PARAMETER 

SECTION 
1 

SECTION 
2 

SECTION 
3 

SECTION 
4 

DIAMETER  (IN) 

WALL  THICKNESS  (IN) 

LENGTH  (IN) 

MOD.  OF  ELAST.  (PSI) 

POISSON'S  RATIO 

WT.  DENSITY  (LB/FT3) 

2.0 
.125 

80.0    , 
30  x  10° 
.30 
490.0 

2.0 
.12  5 

50.0    , 
30  x  10° 
.30 
490.0 

2.0 
.125 

31.6    , 
30  x  10° 
.30 
490.0 

2.0 

.125 

60.0 

30  x 

.30 

490.0 

10' 


DUAL  ANALYSIS  RESULTS* 
MODE 


1 
2 
3 
4 
5 
6 


STARTING  AT 

STARTING  AT 

LEFT  END 

RIGHT  END 

%    DIFFERENCE 

110.5015347 

110.5015347 

.00000 

276.6928874 

276.6928874 

.00000 

314.3597336 

314.3597336 

.00000 

417.0897739 

417.0897739 

.00000 

463.5288123 

463.5288123 

.00000 

675.6526889 

675.6526889 

.00000 

*  SHEAR  DEFLECTION  AND  ROTARY  INERTIA  ARE  CONSIDERED  IN  THIS 
ANALYSIS. 

t  COORDINATES  SHOWN  IN  DIAGRAM  ARE  IN  INCHES. 
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d.   COMPLEX  SYSTEM 

END  CONDITIONS: 

MAIN  MEMBER 
BRANCH 


FIXED-FIXED 
FIXED 


(4>0t2S,-Z5) 


SYSTEM  PARAMETERS 
CALL  SECTIONS) : 

DIAMETER  (IN) 
WALL  THICKNESS  (IN) 
ELASTIC  MODULUS  (PSI) 
POISSON'S  RATIO 


DENSITY    (LB/FT" 


tSJSVY. 


(.0,0,0)    (20,0.0) 


X 


(35.0,0) 


2.25 
.035 
30    x 
.30 
491.5 


10' 


(/&0,20,ZO) 


(60,£G,40) 


(90,-<?tSS) 


(/00,~20,<,o) 


DUAL  ANALYSIS  RESULTS: 


MODE 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 


STARTING  AT 
LEFT  END 

73.0350449 

102.7404907 

142.5065200 

201.7233189 

239.1786608 

327.8210841 

331.8165320 

364.5803428 

428.5622381 

520.1734071 

542.4026997 

565.0762600 

574.0701684 

648.8240390 

754.8640271 

1052.0436077 

1210.9202599 

1250.1288524 


STARTING  AT 
RIGHT  END 

73.7545503 

102.9266533 

146.3168530 

201.3564078 

239.1641100 

321.1763486 

337.3328726 

364.1833110 

429.1841008 

521.2514969 

544.3290021 

564.3679439 

571.0709491 

651.2831922 

759.4399716 

1052.3129250 


%    DIFFERENCE 


1210 
1247 


7321375 
3298693 


985 

181 

674 

182 

006 

,027 

,662 

,109 

.145 

.207 

.355 

.125 

.522 

.379 

.606 

.026 

.016 

.224 
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APPENDIX  G 
SAMPLE  PROBLEM 

1.  Piping  System 

The  piping  used  in  this  example  is  identical  to  the 
complex  system  shown  in  Appendix  F.   It  includes  several  corners 
and  bends  and  has  one  branch  junction  point.   The  boundary 
conditions  are  fixed  for  the  two  ends  of  the  main  member  and 
the  remote  branch  end. 

2.  Data  Cards 

Shown  in  Fig.  G.l   is  the  arrangement  of  data  cards 
as  they  would  appear  following  the  program  deck.   Note  that 
Card  #1  specifies  one  system  to  be  analyzed,  and  Card  #2 
asks  for  four  mode  frequencies  in  the  analysis.   Since  the 
curved  subsection  override  is  not  utilized,  the  program  will 
go  ahead  and  determine  the  number  of  curved  pipe  subdivisions 
according  to  the  length-to-diameter  ratio.   Shear  deflection 
and  rotary  inertia  are  not  included  in  this  problem  analysis. 
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Figure  G.l:   Data  Deck  for  Sample  Program 
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************************************************************************* 

*  * 

*  PROGRAM   VIBREL  CD. RUDOLF  * 

*  * 

************************************************************************* 

PROGRAM  VIBREL  COMPUTES  THE  NATURAL  FREQUENCIES  OF  VIBRATION 

OF  ANY  RANDOMLY  ARRANGED  3-01  ME NS IONAL  PIPING  SYSTEM 
PIPING  HANGERS  AND  PROJECTIVE  COMPLEXITIES  ARE  NOT  INCLUDED 

************************************************************************** 

PROBLEM  #  1 

************************************************************************** 

THE    EFFECTS    OF    SHEAR    DEFLECTION    AND    ROTARY     INERTIA    ARE    NOT    CONSIDERED     IN    THIS    PROBLEM 
************************************************************************** 


NPUT 

DATA: 

X 

Y 

Z 

RAD  OF 

ELASTIC 

POI SSONS 

WALL 

LEFT  BC 

RIGHT  BC 

10 

COORD 

COORD 

COCRO 

CURV 

MODULUS 

RATIO 

0.0. 

THICK 

HT/VOL 

CODE 

CODE 

1 

0.0 

0.0 

0.0 

1 

1 

3 

20.0 

0.0 

0.0 

0.0 

0.300D  08 

.30 

2.25 

0.0350 

491.5 

5 

35.0 

0.0 

0.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

5 

60.0 

25.0 

-25.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

3 

60.0 

25.0 

10.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

4 

60.0 

25.0 

40.0 

3.0 

0.0 

.0 

0.0 

0.0 

0.0 

6 

90.0 

-9.0 

55.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

2 

140.0 

20.0 

-10.0 

1 

4 

140.0 

20.0 

20.0 

3.0 

0.0 

.0 

0.0 

0.0 

0.0 

5 

110.0 

20.0 

20.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

8 

90.0 

-9.0 

55.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

8 

100.0 

-20.0 

60.0 

0.0 

0.0 

.0 

0.0 

0.0 

0.0 

*******  ** ******  ****** $*****************#*«***************************** *** 

PROPERTIES  AND  GEOMETRY: 


SECT 

ELASTIC 

POISSONS 

WALL 

LENGTH  OR 

NR  OF  SUBSEC 

NR 

ID 

MODULUS 

RATIO 

WT/VOL 

O.D. 

THICK 

RAD  OF  CURV 

ALPHA1 

GAMMA 

OR  ALPHA 

1 

1.0 

0.300D 

08 

.30 

491.5 

2.25 

0.0350 

20.000 

0.0 

0.0 

0.0 

2 

3.0 

0.300D 

08 

.30 

491  .5 

2.25 

0.0350 

15.000 

-2.356 

0.955 

0.0 

3 

3.0 

0.300D 

08 

.30 

491.5 

2.25 

0.0350 

43.301 

2.094 

2.186 

0.0 

4 

1.0 

0.300D 

08 

.30 

491.5 

2.25 

0.0350 

35.000 

0.0 

0.0 

0.0 

5 

1.0 

0.300D 

08 

.30 

491  .5 

2.25 

0.0350 

27.833 

0.0 

0.0 

0.0 

6 

2.0 

0.300D 

08 

.30 

491.5 

2.25 

0.0350 

3.000 

1.508 

1.251 

2.000 

7 

5.0 

0.300D 

08 

.30 

491.5 

2.25 

0.0350 

45.592 

0.0 

0.0 

0.0 

8 

1.0 

0.300D 

08 

.30 

491  .5 

2.25 

0.0350 

27.000 

0.0 

0.0 

0.0 

9 

2.0 

0.3000 

08 

.30 

491  .5 

2.25 

0.0350 

3.000 

0.0 

1.571 

2.0J0 

10 

3.0 

0.300D 

08 

.30 

491  .5 

2.25 

0.0350 

27.000 

2.450 

1.156 

0.0 

11 

7.0 

0.300D 

08 

.30 

491  .5 

2.25 

0.0350 

49.659 

-0.372 

-1.177 

-0.  847 

12 

4.0 

0.300D 

08 

.30 

491.5 

2.25 

0.0350 

15.684 

0.0 

0.0 

0.0 

♦♦♦J********************************************************************** 
BOUNDARY  CONDITION  COCES: 

LEFT  B.C.  CO0E=  1 
RIGHT  B.C.  C0DE=1 
BRANCH  #  1  B.C.  C0DE=  1 

************************************************************************** 
PIPING  SYSTEM  NATURAL  FREQUENCIES: 

MODE    FREQUENCYIRAD/SEC)     FREQUENCY  DETERMINANT 

1  73.0350449  0.12278D-22 

2  102.7404907  -.369030-21 

3  142.5065200  0.85913D-21 

4  201.7233189  0.52424D-22 
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GRAPH    OF    THE    VALUE     OF    FREQUENCY    OEIERHIN1NT    VS.     FREQUENCY 
-2.538E-17  3.585E-19 


2.610E-17 


3.8  09E    01 


5.079E    0! 


6.349E     01 


8.888E     01 


1.143E    02 


1.270E    02 


1.651E    02 


1.778E    02 


3.809E     01 


5.079E     01 


6.3-.9E    01 


8.888E     01 


1.143E    02 


1.397E    02 


1.524E    02 


1.905E    02 


2.032E    02 
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